/*
 * GeoReactor.cpp
 *
 *  Created on: 18 Oct 2011
 *      Author: allan
 */

#include "GeoReactor.h"

// GeoReact includes
#include <Geochemistry/GeoSystem.h>

GeoReactor::GeoReactor()
: algorithm("RungeKutta4")
{}

GeoReactor::GeoReactor(const CannonicalSystem& cannonical_system)
: cannonical_system(cannonical_system), algorithm("RungeKutta4")
{
	kinetic_solver.Initialiaze(cannonical_system);
	kinetic_solver.SetControlSolution(0, 1.0E-6);
	kinetic_solver.SetAlgorithm(OdeAlgorithm::RungeKutta4);
}

void
GeoReactor::SetMinimumTimeStep(double dt_min)
{
	kinetic_solver.SetMinimumTimeStep(dt_min);
}

void
GeoReactor::SetMaximumTimeStep(double dt_max)
{
	kinetic_solver.SetMaximumTimeStep(dt_max);
}

void
GeoReactor::SetAlgorithm(const string& algorithm)
{
	OdeAlgorithm ode_algorithm;

	map<string, OdeAlgorithm> ode_algorithms = 
	{
		{"RungeKutta2"                , OdeAlgorithm::RungeKutta2},
		{"RungeKutta4"                , OdeAlgorithm::RungeKutta4},
		{"RungeKuttaFehlberg45"       , OdeAlgorithm::RungeKuttaFehlberg45},
		{"RungeKuttaCashKarp45"       , OdeAlgorithm::RungeKuttaCashKarp45},
		{"RungeKuttaPrinceDormand8"   , OdeAlgorithm::RungeKuttaPrinceDormand8},
		{"RungeKuttaImplicitGaussian1", OdeAlgorithm::RungeKuttaImplicitGaussian1},
		{"RungeKuttaImplicitGaussian2", OdeAlgorithm::RungeKuttaImplicitGaussian2},
		{"RungeKuttaImplicitGaussian4", OdeAlgorithm::RungeKuttaImplicitGaussian4},
		{"BulirschStoerImplicit"      , OdeAlgorithm::BulirschStoerImplicit},
		{"MultistepAdams"             , OdeAlgorithm::MultistepAdams},
		{"MultistepBDF"               , OdeAlgorithm::MultistepBDF}
	};
	
	auto iter_algorithm = ode_algorithms.find(algorithm);;
	
	ode_algorithm = (iter_algorithm != ode_algorithms.end()) ? (iter_algorithm->second) : OdeAlgorithm::RungeKutta4;
	
	kinetic_solver.SetAlgorithm(ode_algorithm);
}

void
GeoReactor::SetAdaptiveTimeStepControl(double eps_abs, double eps_rel)
{
	kinetic_solver.SetControlSolution(eps_abs, eps_rel);
}

const double
GeoReactor::GetMinimumTimeStep() const
{
	return kinetic_solver.GetMinimumTimeStep();
}

const double
GeoReactor::GetMaximumTimeStep() const
{
	return kinetic_solver.GetMaximumTimeStep();
}

const string
GeoReactor::GetAlgorithm() const
{
	return algorithm;
}

const CannonicalSystem&
GeoReactor::GetCannonicalSystem() const
{
	return cannonical_system;
}

void
GeoReactor::TotalCompositionPrimarySpecies(const VectorXd& n, VectorXd& uj)
{
	uj = cannonical_system.TotalCompositionPrimarySpecies(n);
}

void
GeoReactor::TotalEquilibrate(const VectorXd& ujnk, double T, double P, VectorXd& n)
{
	kinetic_solver.TotalEquilibrate(ujnk, T, P, n);
}

void
GeoReactor::ReactStep(KineticStatus& status, double tend, double T, double P, VectorXd& n)
{
	kinetic_solver.ReactStep(status, tend, T, P, n);
}

void
GeoReactor::ReactStep(KineticStatus& status, double tend, GeoSystem& geosystem)
{
	  double T = geosystem.GetTemperature();
	  double P = geosystem.GetPressure();
	VectorXd n = geosystem.GetComposition();
	
	ReactStep(status, tend, T, P, n);
	
	geosystem.SetState(T, P, n);	
}

void
GeoReactor::React(double tbegin, double tend, double T, double P, VectorXd& n)
{
	kinetic_solver.React(tbegin, tend, T, P, n);
}

void
GeoReactor::React(double tbegin, double tend, GeoSystem& geosystem)
{
	  double T = geosystem.GetTemperature();
	  double P = geosystem.GetPressure();
	VectorXd n = geosystem.GetComposition();
	
	React(tbegin, tend, T, P, n);
	
	geosystem.SetState(T, P, n);
}

